Linear response of doped graphene sheets to vector potentials 
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A two-dimensional gas of massless Dirac fermions (MDFs) is a very useful model to describe 
low-energy electrons in monolayer graphene. Because the MDF current operator is directly propor- 
tional to the (sublattice) pseudospin operator, the MDF current-current response function, which 
describes the response to a vector potential, happens to coincide with the pseudospin-pseudospin 
response function. In this work we present analytical results for the wavevector- and frequency- 
dependent longitudinal and transverse pseudospin-pseudospin response functions of noninteracting 
MDFs. The transverse response in the static limit is then used to calculate the noninteracting or- 
bital magnetic susceptibility. These results are a starting point for the construction of approximate 
pseudospin-pseudospin response functions that would take into account electron-electron interac- 
tions (for example at the random-phase-approximation level). They also constitute a very useful 
input for future applications of current-density-functional theory to graphene sheets subjected to 
time- and spatially-varying vector potentials. 

PACS numbers: 71.10.-w,71.45.Gm 
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I. INTRODUCTION 

Graphene, a monolayer of Carbon atoms packed in a 
2D honeycomb lattice, is a recently realized ambipolar 
gapless semiconductor that has attracted enormous in- 
teres1jii2ii2ii^. Electrons in graphene are described at low 
energies by a spin-independent massless Dirac Hamilto- 
nian, which ultimately originates from the non-Bravais 
nature of the 2D honeycomb lattice. The two inequiva- 
lent sites in the unit cell of this lattice are analogous to 
the two spin orientations of a spin- 1/2 particle along the 
-\-z and —z directions (the z axis being perpendicular to 
the graphene plane). This observation opens the way to 
an elegant description of electrons in graphene as parti- 
cles endowed with a pseudospin dcgrce-of-freedo m^i^i^ i^ 
(in addition to the regular spin and valley degrees of 
freedom which play a passive role here). This quantum 
degree-of-freedom has a number of very intriguing impli- 
cations on the electronic properties of this material, most 
of which have been reviewed in Refs. 

Graphene offers a unique example of a new paradigm 
in condensed matter physics: a truly 2D non-Galileian in- 
variant electron hquid (see for example Refs. HHQIIH). 
This non-Galileian invariant nature of graphene is linked 
to non-trivial many-body renormalizations of various 
properties of doped graphene sheets^., such as the plas- 
mon dispersion relation and the optical conductivity^, 
even at very long- wavelengths. Both these properties are 
controlled by the linear response function of 2D MDFs to 
a vector potential, i.e. by the current-current response 
function. Because in the case of MDFs a vector potential 
couples to the pseudospin-density-fluctuation operator, 
this response function happens to coincide, apart from 
a trivial proportionality factor, with the pseudospin- 
pseudospin linear-response function. 

In this work we present analytical results for the 
wavevector- and frequency-dependent longitudinal and 



transverse pseudospin-pseudospin response functions of 
noninteracting 2D MDFs. 

Our paper is organized as follows. In Section [TT] we 
present the model and some basic definitions. In Sec- 
tion |TTT] we present analytical results for the wavevector- 
and frequency-dependent longitudinal and transverse 
pseudospin-pseudospin response functions of a noninter- 
acting undoped 2D MDF system. In Section ITVl we report 
corresponding results for the doped system. In Sect. |V] 
we use the analytical results on the transverse response 
to calculate the orbital magnetic susceptibility of non- 
interacting MDFs. Finally, in Section IVII we present a 
brief summary of our main results and mention their 
usefulness for the construction of (i) response functions 
that would take into account electron-electron interac- 
tions, and (ii) exchange-correlation functionals needed 
in applications of current-density-functional theorjsii to 
graphene sheets in the presence of inhomogeneous vector 
potentials. Finally, two appendices report cumbersome 
technical details and calculations. 



II. THE MODEL AND BASIC DEFINITIONS 

Graphene's honeycomb lattice has two-atoms per unit 
cell and its tt- valence band and tt* -conduction band touch 
at two inequivalent points, K and K' , in the honeycomb 
lattice Brillouin-zone. The energy bands near e.g. the K 
point arc described at low energies by a spin-independent 
massless Dirac Hamiltonian (?i = 1) 



(1) 



k,a,l3 



where "0^ ^ {'i'k,a) creates (destroys) an electron with mo- 
mentum k and sublattice pseudospin a and cr = (cr^, tr^) 
is a vector constructed with two Pauli matrices {cr%i = 



2 



x,y}, which act on the sublattice pseudospin degree-of- 
freedom. Because we are interested in the hnear-response 
functions to smoothly-varying vector potentials, to which 
different spins and valleys contribute independently, in 
Eq. IT]) we have retained only sublattice degrees of free- 
dom. As emphasized repeatedly in the literature (see e.g. 
Refs. 130112), because of the presence of an infinite sea 
of negative-energy states, the MDF model ([1]) must be 
accompanied by an ultraviolet cut-off for the wavevector 
integrals, fc^ax- 

The response to a vector potential A{q, uj) is controlled 
by the current-current linear-response function Xjjili 
which is defined by the usual Kubo producti^ 

^ - lim / di([.4(i),B(0)])e-*e-*.(2) 

The current-density operator jq can be easily found from 
the continuity equation. The density operator pq corre- 
sponding to the Hamiltonian TYd reads 



are two-component pseudospinors. Here A = -1-1 labels 
the conduction band and A = — 1 the valence band and 
Lpk is the angle between k and the x axis, which physi- 
cally denotes the momentum-dependent phase difference 
between wavefunction amplitudes on the A and B sub- 
lattices of graphene's honeycomb lattice. The matrix- 
element factor on the second line of Eq. ^ is given by 



l(XA(fc)|a.|xA'(fc + '7))l' = 



1 + AA' cos((^fe -I- ipk+q) 



Pq =^'Pi-q,a^k,c. ' 



(3) 



(8) 

Note that, because of the continuity equation ([1]), the 
longitudinal pseudospin-pseudospin response function is 
related to the density-density response function Xpp- I^i 
our case this relations reads 

Xpp{^M = ^(['5-q,/5-q]) + --^Xo-a-{qX,U}) . (9) 

This formula holds only at finite w. The first term on 
the r.h.s. of Eq. ([5]), which is an anomalous commutator 
because of the presence of the infinite sea of negative 
energy states, is 'purely real and must be handled with 
great care^ii^. As discussed in Ref. 0, it is easy to show 
that in the noninteracting limit 



and obeys the usual continuity equation 

^^tPq = [Pq,^D] =q jq , 



(4) 



with the current-density operator jq that has a rather 
unusual form2, 



(5) 



The current-density operator for MDFs is proportional 
to the pseudospin-density operator. Due to Eq. ([5]) thus, 
the current-current response function Xjj [Qi '-^) is propor- 
tional to the pseudospin-pseudospin response function. 

In what follows we will calculate the noninteract- 
ing longitudinal pseudospin-pseudospin response func- 
tion, i.e. Xcr^cTj; (^i, with q oriented along the x di- 
rection, q = qx, and the transverse response function, 
Xctx'o-:! ((jy, 1^) with q oriented along the y direction, 
q = qy. These are given byi^ (per spin and per valley) 



— {[^q,P-q] 

vq ^ 



(10) 



where £ma.x/v — fcmax IS an ultraviolet wavevector cut-off. 
Eq. (jlOp is valid to leading order in the limit fcmax ^ oo. 



III. UNDOPED CASE 



We first calculate the response functions x^^^<t^ of the 
undoped system, i.e. the system in which the Fermi en- 
ergy lies at the Dirac point. In this case only the band 
£k,- is full with an occupation factor n^^^_ = 1, while the 
upper band ek.+ is empty (all necessary technical details 
are summarized in Appendix |^ . 

In the longitudinal channel, for what stated at the end 
of Sect.[ni we find the following relation: 
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V q 



^mxf;\q,Lu) (11) 



where 



.(0) 



(q,^) = ^ lim VJ] 



,(0) 



,(0) 



'^fc.A ''^k+q,\' 



X |(xA(fe)k,|xA'(fc + g)>P , 



(6) 



where S is the area of the system, ek,\ — Xvk are single- 
particle Dirac energies, n^"^ are noninteracting band- 
occupation factors and 



^■^^^^ " V2 ( ^6*'^'° 



xf;\l,u.) = -sgn(c.)f-^i^=I=|i (12) 

is the imaginary part of the well-known density-density 
response function of noninteracting MDFs at half fiU- 
ing"'^^. Using the Kramers-Kronig dispersion relations 
one immediately finds the result for the real part 



(7) 5Rexi°l('?i,^) 



2 +4^5ReX^°;H9,^) , (13) 



Attv V q 
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In the trasverse channel, starting again from Eq. ^ 
but with q = qy, we find 

(Oul r ~ \ I ^0(^^^-^'^g^) 

9m XaJMy^"^) = -sgn(w) ^Ju^ - v^q^ . 

Note that the previous equation can also be written as 

v^q 

This relation reminds us of Eq. pT|) . the only differ 



(15) 



ence being that in the transverse case we have the factor 
{up — q^) / {v^ q^) rather than the "gauge-invariance fac- 
tor" Lo^/^v^q^). This simple relation between the trans- 
verse pseudospin-pseudospin response function and the 
density-density response function, however, holds only in 
the undoped case. 

Finally, for the real part of the transverse response 
function, again using the Kramers-Kronig relations, we 
find 



9 2 



5Re xl";'^(9,^) • 
(17) 

In summary, we find that the longitudinal undoped 
pseudospin responses of the model described by ((T|) with 
the ultraviolet cut-off fc,nax are given by 



J 



5Re x^Uqx^uj) 



_£max _ g;^ Q{v^q^-u}^) 
16-D^ - v^q2 



(18) 



while the transverse undoped responses are given by 



?Siex^l%{qy,u;) =- 



Attv 



16t; 



-Qiv^q^-LU^) 



xilt (#,-') = -sgn(^) ^"^' f'^ -eiuj' vV) 



(19) 



16^2 
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Eqs. p8)l and p9)l constitute the first novel results of this 
work. Note that longitudinal and transverse responses 
coincide at g = and that the anomalous terms in the 
first lines of Eqs. (fT8|) and (fTQ]) [first terms on the r.h.s. 
of both equations] are identical. 

Two crucial remarks are in order at this point. Because 
of gauge invariance, the real physical system cannot re- 
spond to a static longitudinal vector potential. More pre- 
cisely, the longitudinal current-current response function 
of the physical system should vanish for cj = and ev- 
ery q, while the transverse one should vanish for cu ~ 
and q ^ 0. After a quick look at Eqs. p8)) and (fTg]) 
one can easily see that this is not what happens to the 
response functions of the model system ([1]), simply be- 
cause of the presence of the cut-off dependent constant 
term — £max/(47rw2). As discussed in Ref. d (see footnote 
33) , this unphysical finite response is due to the fact that 
rigorous gauge-invariance of the model described by Hb 
is broken by the ultraviolet cut-off fcmax- Thus, the static 
response functions of the model system have to be cor- 
rected ad hoc in order to restore gauge invariance: the 
quantity — £max/(47rw2) has to be subtracted away^ from 
the first lines in Eqs. (|18p and On the other hand. 



we would like to remark that in the limit q ~ and 
for finite uj the constant term — £max/(47ru2) is not only 
physical (i.e. it describes the response of the A = — 1 va- 
lence band to a uniform vector potential in the regime in 
which repopulation of states is not allowed^) but also nec- 
essary: when Eq. (fT5)) is substituted inside Eq. this 
term indeed cancels exactly the anomalous commutator 
[by virtue of Eq. (jlOp j. giving a density-density response 
function which is independent of fcmax, as it should^^. 

The pseudospin response functions can be written in a 
more compact form using the following complex functions 



1 



and 



F'Y{q,<^) 

We find 

x(°l(gi(y),c.) = 



9 2 9 



i^L(<Z,c^) 



(20) 



(21) 



47ru2 



z^f^L(T)(g,t^) • (22) 
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FIG. 1: DifTerent regions for the behavior of XctI'itx ('7) '^)- -f^ 
this figure we have introduced dimensionless variables: q — 
q/kp and ui = ui/ep The regions B.l and B.2 are characterized 
by a continuum of interband electron-hole pairs. The regions 

A. 2 and A.3 are characterized by a continuum of intraband 
electron-hole pairs, '^mxala^iq,'^) = in regions A.l and 

B. 3. 



Finally, note that the real part of the long-wavelength 
longitudinal conductivity of the undoped system is 
given by 



3?e crfw) 



? 2 



lini3m (g*, t^) 



e 

16 



(23) 



Section we report, for the sake of completeness, expres- 
sions for x'a}aSqx,uj). 

The contribution to Xo-x''o-x (92;, w) due to doping, 
5xa}a^{qx,uj), is introduced according to 



X!,"J,^ (gi, u;) - (qx, oj) + 5^^,^ (qx, iu) . (24) 

Explicit expressions for the real and imaginary parts of 

Xo-x'ctx (9^1 ^) have been given before in Sect. IIIII Below 
we thus report only expressions for the doping-dependent 

quantity (5xa°k (?*, 

Following Wunsch et alJ^ we introduce the complex 
function 



Gl(z) = z^/z^ - 1 - ln(z + ^z^ - 1) (25) 



and 



UJ± 



2eF ± 
vq 



(26) 



Using G-l{z) and the function Fi^{q,uj) introduced above 
in Eq. ^ we find 



Sx^^Liqx,^) = -|^2^+FL(<z,c.){GL(a.+ ) 

- e(cc;_-l)[GL(co_)-z7r] 

- e(l-^_)GL(-a>_)} . (27) 

We now provide more explicit expressions for 
Sx'a^a^iqx.uj) in terms of real functions. With ref- 
erence to Fig. [1] we find the following analytical results. 



Restoring h and introducing the gsgv = 4 spin-valley 
degeneracy we find the usual "universal" frequency- 
independent valu&i^ SRe cr{uj) = e^/(4?i). 



IV. DOPED CASE 



We now pass to calculate the response functions Xo 
of the doped system, i.e. the system in which the Fermi 
energy lies, for example, above the Dirac point. In this 
case the band e^,- is completely full with the usual occu- 
pation factor v^^_ = 1, while the upper band £u.+ is filled 

only up to the Fermi energy ep = vfcp, where /cp = \f^irn. 
In the case of finite doping, we find more convenient to 
evaluate Eq. ([5]) on the imaginary-frequency axis (i.e. by 
letting w iuj) and then to perform a standard analyti- 
cal continuation to the real-frequency axis (more details 
are given in Appendix [B|) . 



1. Region A.l 



For uj < vq — 2£p: 



^eSx^^],(qx,iu) = - 



UJ £f 

v^q^ 27ru^ 
, ,2 



167ru^ \J v'^q^ — 



aL(q,t^) 



and 3m 5x^a}cj^{qx,uj) = 0. Here 

ai,{q,uj) = arcsin (w+) -I- w+ y/ 1 — ui 



+ arcsin (ti>_ ) + a / 1 — 



qion A.2 



(28) 



A. Longitudinal channel 

Even though, as discussed before, the longitudinal re- 
sponse function is determined by the density-density re- 
sponse functio n^^i^^i^^ via the continuity equation, in this 



For Lu < vq and to > \2e-p — vq\: 
3te 6x^J>],Jqx,L.) 



Wv^y/v^q^ -u;^ v^q^ 27rw2 
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and 



5. Region B.2 



3m <5xS!kM,^) = — — ^====CL(g,a>) . (30) 



Here 



bi,{q,uj) — arccos(a;_) — ti>_^l — cj^ (31) 



and 



CL(g, w) =ln[u;+ + Juj\ - 1 - uj+ JujX - 1 . (32) 



3. Region A. 3 



For u) < vq and lo < 2eF — vq: 
and 



:dL(g,Cc) . (34) 



Here 



^^1+ + wcj^ - 1 

dL(g,t^) = In ( I - tJ+A/cji - 1 

w_ + i/w^ — 1 



+ w_-\/a;l — 1 



(35) 



4- Region B.l 



For u; > 2ep + wg: 



, ,2 



peL(g,w) 



(36) 



and 3m (5x1° (gi, w) = 0. Here 



eL(g,a;) = In 



w I + a/o;! - 1 



Lj_ — 1 — 



- w+a/cj^ - 1 



UJ-\ LU_ — 1 



(37) 



For LO > vq, lo > 2eF — vq, and lo < 2eF + vq: 



Sx^^},Jqx,Lo) = - 



Lo^ ep 
v'^q^ 2t:v^ 



IGttv'^ ^ LlP- ~ v'^q' 



(38) 



and 



Sm 5xi°L(9*,t^) 



IGtTW^ y^LO^~V^' 



F.9L(g,w) 



(39) 



Here 



(33) /L(g,a;) = In + Jc^^ _ ^ _ cc>+a/cj2 _ ^ (40) 



and 



ghiq,Lo) = arccos (lo^ ) — uj-\J 1 — u"^ . (41) 



6. Region B.3 



For LO > vq and w < 2eF — vq: 



v'^q^ 2'KV^ 



IGttu^ ^ LoP- ~ v'^q^ 



hi,{q,uo) 



and 



"^m 6x''°}„Jqx,uj) 



IGv'^y^Lo'^ — v'^q^ 



(42) 



(43) 



Here 



hi^{q,uj) = In 



L0+ + \l Uj\ - 1 



UJ- + X Uji - 1 



aj+A/cji - 1 



U!-\ L0_ — 1 



(44) 



We would now like to remind the reader that in 
the region B.3 of the {q,io) plane the interacting sys- 
tem possesses a collective plasmon mode LOpi = LOpi{q), 
which, within the simple random phase approximation 
(RPA)i^iiSiiiii^, can be found by solving the equation 



l-v,^exi°J{q,io)=0 , 



(45) 
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where Vq = 27re^/(eg) is the 2D Fourier transform of In what follows, we report only expressions for 



the Coulomb potential (e being an average dielectric con- 
stant that depends on the environment surrounding the 
graphene flake). Using Eq. ^ we can write Eq. (|^5)) in 
the more appealing form 



1 



V q 







(46) 



Using the microscopic expression for SRe x^^^^^nil^^^) 
Eq. (|18p and Eq. (|^^ , expanding the expression in square 
brackets in Eq. (|46p in powers of q/k-p up to fourth order, 
introducing the gsg^ = 4 spin-valley degeneracy factor, 
and restoring for a moment h, we find that the plasmon 
dispersion is given by 



(47) 

where we have introduced the density-dependent "plas- 
mon mass" 



4?™ 



(48) 



the fine structure constant acc = e^/(^'^e); and the 
Thomas- Fermi screening vector fcxF = ffsffvOeefep- 

Note that mpi tends to the bare electron mass in vac- 
uum m if we use the parabolic energy-momentum dis- 
persion s-p = h^kp/{2m) (while still using gsgv = 4). 
In graphene, however, rupi = kkp/v, which (i) is cx y/n 
and (ii) trivially depends on Planck's constant h because 
the Fermi energy in this material depends linearly on h 
(rather than quadratically) . As explained in Ref . [^, be- 
cause the MDF model Hamiltonian is not invariant under 
an ordinary Galileian boost, RPA is not exact for inter- 
acting systems of MDFs even in the limit q ^ (con- 
trary to what happens in the conventional 2D parabolic- 
band electron gas^^, where for q the plasmon is pro- 
tected from many-body renormalizations by Galileian in- 
variancc). When interactions between MDFs are treated 
beyond RPA the plasmon mass acquires a non-trivial 
density and coupling-constant dependence^. 

Last but not least, note that, at odds with what stated 
in Ref. 17, the first subleading correction to the RPA 
plasmon dispersion [second term in round brackets in 
Eq. (|T7l) ] can change sign: it is positive (like in the con- 
ventional 2D parabolic-band electron gas^'^) at weak cou- 
pling, but becomes negative for a^c > \/3/2. 



B. Transverse channel 



In this Section we report explicit expressions for 

the transverse response Xo^-^i^xilVi^)- As done in 
Sect. IIV Al we introduce the contribution due to doping, 

<5x^Sx('7y,t^), according to 

xSt^ {qy, u;) = xtl ^) + 6x^°L in, CO) . (49) 



'5x^S.('7y,'^)- 

Similarly to what done in Sect. IIV"S] we introduce the 
complex function 



Gt(z) = zy/z^ - 1 + ln(z -|- y/ - 1) . (50) 

Using Gt{z) and the function FT^{q,Lo) introduced above 
in Eq. ^ we find 



v'^q^ 2Trv^ 



FT(q,cj){GT(t^+) 



- e(w_ - 1)[Gt(w_) +i7r] 

- e{l-uj-)GT{-uJ-)} . (51) 

Once again, we provide below more explicit expressions 
for 6xlj](7^{qy,i^) in terms of real functions. 



A.l 



For Lu < vq — 2eF: 



£f y/v'^q'^ - 



(52) 

and 3m (Jxa^CT^ (f/y, t^) = 0. Here 

aT{q,uj) — arcsin (ti^+) — lu+ \J ^ ^ ^+ 



-|- arcsin (cj_) — W-Wl — w 



2. Region A. 2 



For Lu < vq and to > \2e-F — vq\: 



5x^^},(qy,Lo) = - 



Il2g2 2TrV^ 



y/v^Cp 



Wttv^ 



-bT{q,io) (53) 



and 



q — UJ 



Sx^^l^{qy,co) = - ' CT(g,^) ■ (54) 



Here 



6t (q,!^) — arccos((jj_ ) + ui-yjl — uPi (55) 



and 



ct((J, w) = In CJ+ + v /o;^ - 1 -^uj+Juj\-\ . (56) 
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3. Region A. 3 
For Lo < vq and w < 2eF — vq: 



^e5x^Sl{qV.uj) = -^L^ + -—-1-^ (57) 



16^2 



and 



5m5x^°L(#,c.) = -V^-^dT(g,c^) . (58) 



Here 



(ix('?j'^) — In 



Ijj I + . /cj|l — 1 



CiJ_ + \/ UJ1_ — 1 



UJ-\ UJ_ — 1 



(59) 



4- Region B.l 



For uj > 2eF + wq: 
and 3m <5xo-x''o-x (9^) '^) = 0. Here 



^^^^ — 2 T7 — 2 eT(g,w) 



(60) 



eT(g,w) = In I ; I +U!+Juj^ - 1 



LO^\luj_ — 1 . 

5. Region B.2 
For LO > vq, uj > 2eF ~ vq, and uj < 2£f + vq- 



(61) 



and 



(62) 



\/ up — v'^q^ \J up ~ v'^q^ 



low lOTTW^ 



Here 



(63) 



/t(q,w) = In UJ+ + Jul ~ 1 +u;+Ju;l ~ 1 (64) 



and 



(7t(q, w) = arccos ((jj_) + w^y^l — . (65) 





2 3 



4 5 6 



FIG. 2: The imaginary part of the longitudinal pseudospin 
response function, — SJm x12o-^(Q' ~ 92;, tj) [in units of the 
density-of-states at the Fermi level, /^(O) = ep/(27ru^)], as a 
function of tj/ep. a) g = 0.5 fcp, b) g = 1.0 fcp, c) g = 2.0 fcp, 
and d) g = 3.0 ftp. Singularities at a; = wg are clearly visible^^. 



6. Region B.3 
For UJ > vq and u; < 2eF — vq: 



and 



^JuP--'\Pq^ 
16f 2 



(66) 
(67) 



Here 



hT{q,Uj) = 




UJ-\ UJ 



(68) 



In Figs. [Un] we have reported plots of the longitudinal 
and transverse response functions. Note that the trans- 
verse response function goes always to zero at w = vq: 
in this case, indeed, the angle between fc and q is either 
zero (intraband term) or tt (interband term). Recalling 
that in the transverse case q is directed along y, we thus 
find that either (pfe+q — ipk = /'i (intraband term) or 
Lpk+q = — </?fc — 7i'/2 (interband term): in both cases the 
matrix element in Eq. ([5]) vanishes. Finally, note that 

both 5to Xo-x'ct^ (9^) ^) a-iid 5m x^o}a^ {qy, w) diverge lin- 
early for UJ 3> vq. 
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uj/sf uj/sp ^/ef w/sf 



FIG. 3: The imaginary part of the transverse pseudospin re- 
sponse function, —9m x12ct^(9 ~ IVj'^) [i^ units of !^(0)], as 
a function of a;/eF- a) g = 0.5 fcp, b) g = 1.0 fcp, c) g = 2.0 ftp, 
and d) g = 3.0 fcp. 



FIG. 5: The real part of the transverse pseudospin response 
function, — Ke xi2<T="('? ~ IVj^) [i^ units of 1^(0)], as a func- 
tion of ui/ep. a.) q — 0.5 fcp, b) g = 1.0 fcp, c) g = 2.0 fcp, and 
d) g = 3.0 fcp. 
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FIG. 4: The real part of the longitudinal pseudospin response 
function, — Ke XCT°io-^('7 = 'lS;,u}) [in units of i^(0)], as a func- 
tion of uj/ef. a) g = 0.5 fcp, b) g = 1.0 fcp, c) g = 2.0 fcp, and 
d) g = 3.0 fcp. The filled circles in panels a) and b) mean that 

the function — Ke X^^^a^il ~ 1^^^) takes exactly the value 1 
[in units of v{0)] at ui = vq. 



calculated from 



(0) 
-^orb 



y2 2 (0) 

— lim 



(9y,0) 



(69) 



As stressed earlier in Sect. IIIIl before taking the limit in 
Eq. ((69|) the undoped contribution to the static response 



xi^o-ic (gy, 0) has to be regularized to restore gauge in- 
variance by subtracting the cut-off dependent constant 
term — £max/(4vru^). A simple inspection of the equa- 
tions reported in the previous section allows us to write 
a compact expression for the static transverse response 
(as usual, per spin and valley): 



Act a. 



(qy, 0) = 1 1 - - arcsin [£{q)] 

IbV I TT 



where 



„, . 1 / 2fcF 



2 fcp 

q 



e 1 



2/cp 



1 - 



2 fcp 

q 



(70) 



(71) 



V. DIAMAGNETIC SUSCEPTIBILITY 

The transverse pseudospin response function allows us 
to calculate the orbital magnetization induced in doped 
graphene sheets by a static magnetic field. As discussed 
in detail in Sects. 3.4.3. and 4.5 of Ref. [3, the non- 



interacting orbital magnetic susceptibility Xorb '^^n be Xo-Jcr^{qy,0)/q is proportional to (S(£f). Indeed, if ip{e) 



We thus immediately see that for ep > and for all 
wavevectors q < 2k^, £{q) = 1 and Xo-^^ (^y, 0) = 0. 
This implies that the orbital magnetic susceptibility of 
noninteracting MDFs is zero. However, if ep = 

x'a^a^ {qy, 0) (X g and thus the orbital magnetic suscep- 
tibility diverges in the undoped limit. More precisely, it 
is possible to show that in the limit g — s- the function 

(0) 



9 



is a test function (here e is a shorthand notation for ep), 



q 



PVq/2 ^ 

2 lim / de -— 



1 2 . f2e 

1 arcsm — 

TT \vq 



2 2e 
TT vq 



vie). 



(72) 



Introducing the dimensionless variable x = 2e/ (vq) we 
find 



g^O I 16 Stt 



/ dx arcsin(a;) 

^0 



dx x\f\ — xP- 



>^(e) = — (^(0) , (73) 
ovr 



where e S [0,f(7/2]. In summary, introducing the ^sffv = 
4 degeneracy factor, we find^^ 



(0) 
•^orb 



67r 



-5(£f) 



(74) 



in perfect agreement with Ref . [20- As explained by Mc- 
Clure^, the origin of the large (infinite at T = 0) dia- 
magnetism in undoped graphene can be understood qual- 
itatively. When the magnetic field B is turned on, group 
of states, which were originally distributed in energy, co- 
alesce into Landau levels. In undoped graphene, states 
which had negative energy coalesce into the n = Lan- 
dau level, thus increasing the energy of the system, which 
will respond to the field with a negative orbital suscep- 
tibility. The total energy (per unit area) in the absence 
of the field of the group of electrons which condense into 
the n = Landau level is 



En 



de ev{e)f{e) 



deey{e)[f{e)~f{~e)], (75) 







where Wc = a/2w/£b c>c ^/B is the MDF cyclotron fre- 
quency [£_B = \/c/[eB) being the magnetic length], 
v{e) = |e|/(27rz;^) is the density-of-states at B = 0, and 
/(e) is the Fermi-Dirac distribution function. The ex- 
treme of integrations in the first line of Eq. (j75p ensure 
that the number of states 



+ A 



de v{e) 



A2 



(76) 



can be accommodated into the n = Landau level, which 
has degeneracy eB/{2'Kc) per unit area, i.e. A = lOc/V2. 
In the zero-temperature limit 

c/\/2 

Eq = -TT^ I de e^ 



2ttv^ 



1 



12^/27! 



(77) 



The total energy in the presence of the field is zero 
(because turning on the B field all the electrons con- 
dense into the n — Landau level, which is at zero 
energy). Thus the change in energy with the field is 
AE = -Eo oc B'-^/^. The susceptibility is thus 



(0) 
Aorb 



1 d{AE) 
hm — — -— — cx 

B^a B dB 



^ ' (78) 



which diverges for B ^ Q. 

Introducing the usual Bohr magneton, /iB = eh/ (2mc), 
and restoring ?i, the orbital magnetic susceptibility ([7^ 
can be written as 



(0) 
Aorb 



V 2 / 



37r?i 



<5(£f) 



(79) 



g ^ 2 being the bare electron g-factor in vacuum. We 
recall that the Pauli spin susceptibility of 2D MDFs is 



(0) 



2iTh^v^ 



(80) 



In the conventional 2D parabolic-band electron gas 
xi°b = -XpVS- For 2D MDFs we see that the den- 
sity dependence of x^oX, completely different from that 

( (0) 

of Xp . 



VI. SUMMARY AND CONCLUSIONS 

In summary, we have calculated analytically the lon- 
gitudinal and transverse pseudospin-pseudospin linear- 
response functions of noninteracting massless Dirac 
fermions. As expected, because of the continuity equa- 
tion that relates the density operator with the longitu- 
dinal current operator, the longitudinal response func- 
tion is determined by the density-density response func- 
tion (apart from an anomalous commutator term, which 
has to be handled carefully). We have used the trans- 
verse pseudospin response function to calculate the or- 
bital magnetization induced in doped graphene sheets by 
a static magnetic field, finding perfect agreement with 
earlier calculations'^" based on the explicit use of the 
Landau level structure of two-dimensional massless Dirac 
fermions. 

The results presented in this work constitute a very 
useful starting point for the construction of approximate 
response functions that would take into account electron- 
electron interactions. For example, within the random 
phase approximation^ the pseudospin-pseudospin re- 
sponse functions of the interacting system can be written 
as 



x<T-<T-(q,w) 



X, 



1 



(81) 



Our results are also very useful in view of future appli- 
cations of current-density-functional theoryii to doped 
graphene sheets in the presence of time- and spatially- 
varying vector potentials (see also comments in Sect. II 
of Ref.^). 
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It is easy to see that z+Z- = 1 and that jz+j < 1, which 
implies |z_| > 1. Calculating the residue in (which is 
a third-order pole) we finally find Eq. (jisp . 



2. Real part 



APPENDIX A: DETAILS ON THE ANALYTICAL 
CALCULATION OF THE TRANSVERSE 
UNDOPED RESPONSE FUNCTION 

1. Imaginary part 

In the undoped limit and for uj > 0, the imaginary part 
of the transverse response function is given by 



(2^)2 
1 - C0s(v3fc+q + ipk) 



S {uj — vk — v\k + q\) 
(Al) 



This expression can be easily obtained from Eq. ^ by 
retaining only the interband term corresponding to A = 
+ and A' = — (because uj > 0). The cosine term in the 
second line of Eq. (jAl|) can be easily manipulated to give 



cos(^fc+q + Lpk) = \kTq\ ' ^ ^ 

For a given value of uj the delta function in Eq. (jAl[) gives 
a non-zero contribution to the fc-integration if and only 
if 



UJ — vk = v\k + q\ 
which can be solved for k yielding 



9 ? ? 

UJ^ — V q 



2v{u! + vqsin ipk) 
Performing the integration over k we are left with 



(A3) 



(A4) 



9 9 2 



2tt 



Q{ijj — vq) 
{ivsinipk + vq)'^ 



{uj + vqsin ipk)^ 

(A5) 

The angular integration can be easily performed in the 
complex plane z = exp{iipk). We thus find 

10 {uj/vq + sin ipkf Jc [(z - z+)(z - z_)]3 ' 

(A6) 

where C is the unit-radius circle in the complex plane and 



Zl,2 



Z± 



-ivq 



2„2 



-oj ± ^ LlP- — v^q^ 
vq 



(A7) 
(A8) 



We now pass to calculate JRe Xo-^ctL (w, by using the 
Kramers- Kronig dispersion relation, i. e. 

» (On) . ^ . '^^ r ^'^^xttiqA.u;') 



1 



lim V 



cW^^H^Z. (A9) 



vq 



Performing the change of variable — uj''^ — v'^q^ and 
defining i^ax = V^max - v'^Q^ we find 



8ttv^ 



X V 



Snv^ 
dt 



,(A10) 



where the first term on the r.h.s. diverges in the limit 
(^max ~* oo. Considering that Wmax is a maximum exci- 
tation energy [and thus it is twice the cut-off Emax intro- 
duced after Eq. (fT^l) ]. this term can be re- written as 



8nv^ 



1 - 



vq 



2e„ 



(All) 



In the limit Emax ~* oo this terms thus tends to 
— e,nax/(47rv^). The second integral gives 



V 



dt 



n Q(vq — uj) 

i2_(^2_„2^2y - 2^^!^f=^ 



(A12) 



Summing together these two contributions we find im- 
mediately the final result ^ 



APPENDIX B: DETAILS ON THE ANALYTICAL 
CALCULATION OF THE TRANSVERSE DOPED 
RESPONSE FUNCTION 

In this Appendix we report some details on the analyti- 
cal calculation of the transverse doped response function. 
In the doped case it turns out more convenient to per- 
form the calculation on the imaginary frequency axis, i.e. 
LU iuj in Eq. ([6]): 



(0) _ (0) 



S iuj + X'vk - \v\k + q\ 
1 -I- AA' cos [ipk+q + (ph) 



(Bl) 
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Subtracting the undoped contribution and doing simple 
algebraic manipulations we find 



Sxtl. iiy^ ^) = Jil^ ^) + J* (9: ^) : (B2) 



where 



k-p 



"^(*'^)^^y ^'^^ \ dVkf{k,q,uj) (B3) 



with 



f{k,q,uj) = 



2vk sin^ ifk + vq sin ip^ — {i'^ + "^vk) 



2v'^kqsmipk + lu'^ + v^tp — 2iujvk 

(B4) 

Once again the angular integration can be easily per- 
formed in the complex plane z = exp {i(fk)'- 



1 /•^'^ 

^^1^^) ^ ^ j d,ipkf{k,q,Lj) 



dz 



P{z) 



St:^ v^kq J (2 z'^{z — Zj^){z — z^) ' 



(B5) 



with P{z) = vk{z^ - ly + ivqz{z^ - 1) + 2(iw + 2vk)z^, 



z± = 



-i[v'^q^ — {iujY ~ 2iujvk] 
2v^kq 



It is possible to show that z+z- = —1 and \z^\ < 1. 
Calculating the residues in z = (a second-order pole) 
and in 2 = z_)_ (a first-order pole) and performing the 
integration over k in Eq. (|B3|) we finally find 



Jiq,uj) 



167rt>2 



arcsin(ri) + i^\J 0^ — 1 

(B7) 



where $7 = (2e-p^iu^)l(vq). Using this result in Eq. (|B2p 
and performing a standard procedure of analytical con- 
tinuation one finds the results in Sect. IIVI 
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